#define MAX_SIZE 1111
#define MAX_THREADS 8

*****************************************************************************
*
*	Array generation
*
*****************************************************************************
      subroutine sgen( a, n)
      real a(*)
      integer n
      integer i,ir

      ir = mod(1325 * (irand()+1), 65536)
      ratio = 1.0 / 32768.0
	do i = 1, n
	    ir = mod(ir*3125,65536)
	    a(i) = float( ir ) * ratio - 1.0
	end do
      return
      end

      subroutine dgen( a, n)
      double precision a(*)
      integer n
      integer i,ir

      ir = mod(1325 * (irand()+1), 65536)
      ratio = 1.0 / 32768.0
	do i = 1, n
	    ir = mod(ir*3125,65536)
	    a(i) = float( ir ) * ratio - 1.0
	end do
      return
      end

      subroutine cgen( a, n)
      complex a(*)
      integer n
      integer i,ir
      real re, im, ratio

      ir = mod(1325 * (irand()+1), 65536)
      ratio = 1.0 / 32768.0
	do i = 1, n
	    ir = mod(ir*3125,65536)
	    re = float( ir ) * ratio - 1.0
	    ir = mod(ir*3125,65536)
	    im = float( ir ) * ratio - 1.0
	    a(i) = cmplx( re, im)
	end do
      return
      end

      subroutine zgen( a, n)
      double complex a(*)
      integer ld, n
      integer i,j,ir
      double precision re, im, ratio

      ir = mod(1325 * (irand()+1), 65536)
      ratio = 1.0 / 32768.0
	do i = 1, n
	    ir = mod(ir*3125,65536)
	    re = float( ir ) * ratio - 1.0
	    ir = mod(ir*3125,65536)
	    im = float( ir ) * ratio - 1.0
	    a(i) = cmplx( re, im)
	end do
    
      return
      end

*****************************************************************************
*
*	Array Sum
*
*****************************************************************************
      real function ssum( n, x, incx)
      integer n, incx
      real x(*)
	ssum = 0.0d0
	ix = 1
	do i = 1, n
	    ssum = ssum + x(ix)
	    ix = ix + incx
	enddo
      return
      end

      double precision function dsum( n, x, incx)
      integer n, incx
      double precision x(*)
	dsum = 0.0d0
	ix = 1
	do i = 1, n
	    dsum = dsum + x(ix)
	    ix = ix + incx
	enddo
      return
      end

*****************************************************************************
*
*	1D FFT
*
*****************************************************************************
      subroutine sft1d( plusminus, n, a, inc)
      real a(*)
      real w(MAX_SIZE)
      integer n, plusminus, inc
      
      real u, v
      TPI = 6.28318530717959

      if ( plusminus .eq. -1) then
	w(1) = ssum( n, a, inc)
	l = (n+1)/2
	do k = 2, l
	    u = 0.0
	    v = 0.0
	    do i = 1, n
		u = u + a(inc*(i-1)+1) * cos((k-1) * (i-1) * tpi / float(n))
		v = v - a(inc*(i-1)+1) * sin((k-1) * (i-1) * tpi / float(n))
	    end do
	    w(2*k-2) = u
	    w(2*k-1) = v
	end do
	if( mod(n,2) .eq. 0) then
	    u = 0.0
	    do k = 1, l
		u = u + a(inc*(2*k-1-1)+1) - a(inc*(2*k-1)+1)
	    end do
	    w(n) = u
	end if
	do i = 1, n
	    a(inc*(i-1)+1) = w(i)
	end do
      elseif (plusminus .eq. 1) then
	Print *, 'Sorry NOT Implemented'
      else
	print *, 'Error first argument of SFT1D should -1 or +1'
	stop
      endif

      return
      end

      subroutine dft1d( plusminus, n, a, inc)
      double precision a(*)
      double precision w(MAX_SIZE)
      double precision dsum
      integer n, plusminus, inc
      
      double precision u, v, tpi
      TPI = 2 * 3.14159265358979323846

      if ( plusminus .eq. -1) then
	w(1) = dsum( n, a, inc)
	l = (n+1)/2
	do k = 2, l
	    u = 0.0
	    v = 0.0
	    do i = 1, n
		u = u + a(inc*(i-1)+1) * cos((k-1) * (i-1) * tpi / float(n))
		v = v - a(inc*(i-1)+1) * sin((k-1) * (i-1) * tpi / float(n))
	    end do
	    w(2*k-2) = u
	    w(2*k-1) = v
	end do
	if( mod(n,2) .eq. 0) then
	    u = 0.0
	    do k = 1, l
		u = u + a(inc*(2*k-1-1)+1) - a(inc*(2*k-1)+1)
	    end do
	    w(n) = u
	end if
	do i = 1, n
	    a(inc*(i-1)+1) = w(i)
	end do
      elseif (plusminus .eq. 1) then
	Print *, 'Sorry NOT Implemented'
      else
	print *, 'Error first argument of SFT1D should -1 or +1'
	stop
      endif

      return
      end


*****************************************************************************
*									    *
*	2D FFT								    *
*									    *
*****************************************************************************
      subroutine sft2d( plusminus, n1, n2, w, ldw)
      real w(ldw,n2)
      real save( 2*MAX_SIZE), cp(2*MAX_SIZE,MAX_THREADS)
      integer n1, n2, lda, ldw, plusminus
      

      if (plusminus .eq. -1) then
	call sfft1di( n1, save)
c$doacross local(j), share(cp)
	do j = 1, n2
	    call sfft1d( plusminus, n1, w(1,j), 1, save)
	end do

	call sfft1di( n2, save)
	call sfft1d( plusminus, n2, w(1,1), ldw, save)
	if( mod(n1,2) .eq. 0) then
	    call sfft1d( plusminus, n2, w(n1,1), ldw, save)
	end if

	call cfft1di( n2, save)
c$doacross local(i,j,my)
	do i = 2, n1-1, 2
	    my = 1
c$	    my = mp_my_threadnum()+1
	    do j = 1, n2
		cp(2*j-1,my) = w(i+0,j)
		cp(2*j-0,my) = w(i+1,j)
	    end do
	    call cfft1d( plusminus, n2, cp(1,my), 1, save)
	    do j = 1, n2
		w(i+0,j) = cp(2*j-1,my) 
		w(i+1,j) = cp(2*j-0,my)
	    end do
	end do
      else
	Print *, 'Sorry NOT Implemented'
      end if
      return
      end
      

      subroutine dft2d( plusminus, n1, n2, w, ldw)
      double precision w(ldw,n2)
      double precision save( 2*MAX_SIZE), cp(2*MAX_SIZE,MAX_THREADS)
      integer n1, n2, lda, ldw, plusminus
      

      if (plusminus .eq. -1) then
	call dfft1di( n1, save)
c$doacross local(j)
	do j = 1, n2
	    call dfft1d( plusminus, n1, w(1,j), 1, save)
	end do

	call dfft1di( n2, save)
	call dfft1d( plusminus, n2, w(1,1), ldw, save)
	if( mod(n1,2) .eq. 0) then
	    call dfft1d( plusminus, n2, w(n1,1), ldw, save)
	end if

	call zfft1di( n2, save)
c$doacross local(i,j,my), share(cp)
	do i = 2, n1-1, 2
	    my = 1
c$	    my = mp_my_threadnum()+1
	    do j = 1, n2
		cp(2*j-1,my) = w(i+0,j)
		cp(2*j-0,my) = w(i+1,j)
	    end do
	    call zfft1d( plusminus, n2, cp(1,my), 1, save)
	    do j = 1, n2
		w(i+0,j) = cp(2*j-1,my) 
		w(i+1,j) = cp(2*j-0,my)
	    end do
	end do
      else
	Print *, 'Sorry NOT Implemented'
      end if
      return
      end
      

*****************************************************************************
*
*	3D FFT
*
*****************************************************************************
      subroutine sft3d(plusminus,n1,n2,n3,w,ld1,ld2)
      real w(ld1,ld2,n3)
      real save(2*MAX_SIZE)
      integer n1,n2,n3,ld1,ld2,plusminus

      if( Plusminus .ne. -1) then
	Print *, 'Sorry NOT Implemented'
	return
      endif
c
c   transform along X  first ...
c
      call sfft1di( n1, save)
c$doacross local(i,j,k)
      do k = 1, n3
       do j = 1, n2
	    call sfft1d( plusminus, n1, w(1,j,k), 1, save)
       end do
      end do
c
c   transform along Y then ...
c
      call sfft1di( n2, save)
c$doacross local(i,j,k)
	do k = 1,n3
	    call sfft1d( plusminus, n2, w(1,1,k), ld1, save)
	end do
      if (mod(n1,2) .eq. 0) then
c$doacross local(i,j,k)
	do k = 1,n3
	    call sfft1d( plusminus, n2, w(n1,1,k), ld1, save)
	end do
      end if

      call cfft1di( n2, save)
c$doacross local(i,k)
      do k = 1,n3
	do i = 2, n1-1, 2
	    call csfft1d( plusminus, n2, w(i,1,k), 1, ld1, save)
	end do
      end do
c
c   transform along Z finally ...
c
      call sfft1di( n3, save)
	call sfft1d( plusminus, n3, w(1,1,1), ld1*ld2, save)
	if ( mod(n2,2) .eq. 0) then
	    call sfft1d( plusminus, n3, w(1,n2,1), ld1*ld2, save)
	end if
      if (mod(n1,2) .eq. 0) then
	    call sfft1d( plusminus, n3, w(n1,1,1), ld1*ld2, save)
	if ( mod(n2,2) .eq. 0) then
	    call sfft1d( plusminus, n3, w(n1,n2,1), ld1*ld2, save)
	end if
      end if

      call cfft1di( n3, save)
c$doacross local(i,j,k)
	do j = 2, n2-1, 2
		call csfft1d( plusminus, n3, w(1,j,1), ld1, ld1*ld2, save)
	    if (mod(n1,2) .eq. 0) then
		call csfft1d( plusminus, n3, w(n1,j,1), ld1, ld1*ld2, save)
	    end if
	end do
c$doacross local(i,j)
      do j = 1,n2
	do i = 2, n1-1, 2
	    call csfft1d( plusminus, n3, w(i,j,1), 1, ld1*ld2, save)
	end do
      end do

      return
      end

      subroutine dft3d(plusminus,n1,n2,n3,w,ld1,ld2)
      double precision w(ld1,ld2,n3)
      double precision save(2*MAX_SIZE)
      integer n1,n2,n3,ld1,ld2,plusminus

      if( Plusminus .ne. -1) then
	Print *, 'Sorry NOT Implemented'
	return
      endif
c
c   transform along X  first ...
c
      call dfft1di( n1, save)
c$doacross local(i,j,k)
      do k = 1, n3
       do j = 1, n2
	    call dfft1d( plusminus, n1, w(1,j,k), 1, save)
       end do
      end do
c
c   transform along Y then ...
c
      call dfft1di( n2, save)
c$doacross local(i,j,k)
	do k = 1,n3
	    call dfft1d( plusminus, n2, w(1,1,k), ld1, save)
	end do
      if (mod(n1,2) .eq. 0) then
c$doacross local(i,j,k)
	do k = 1,n3
	    call dfft1d( plusminus, n2, w(n1,1,k), ld1, save)
	end do
      end if

      call zfft1di( n2, save)
c$doacross local(i,k)
      do k = 1,n3
	do i = 2, n1-1, 2
	    call zdfft1d( plusminus, n2, w(i,1,k), 1, ld1, save)
	end do
      end do
c
c   transform along Z finally ...
c
      call dfft1di( n3, save)
	call dfft1d( plusminus, n3, w(1,1,1), ld1*ld2, save)
	if ( mod(n2,2) .eq. 0) then
	    call dfft1d( plusminus, n3, w(1,n2,1), ld1*ld2, save)
	end if
      if (mod(n1,2) .eq. 0) then
	    call dfft1d( plusminus, n3, w(n1,1,1), ld1*ld2, save)
	if ( mod(n2,2) .eq. 0) then
	    call dfft1d( plusminus, n3, w(n1,n2,1), ld1*ld2, save)
	end if
      end if

      call zfft1di( n3, save)
c$doacross local(i,j,k)
	do j = 2, n2-1, 2
		call zdfft1d( plusminus, n3, w(1,j,1), ld1, ld1*ld2, save)
	    if (mod(n1,2) .eq. 0) then
		call zdfft1d( plusminus, n3, w(n1,j,1), ld1, ld1*ld2, save)
	    end if
	end do
c$doacross local(i,j)
      do j = 1,n2
	do i = 2, n1-1, 2
	    call zdfft1d( plusminus, n3, w(i,j,1), 1, ld1*ld2, save)
	end do
      end do

      return
      end

